{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Ничипоренко Александр Владимирович\n",
    "    \n",
    "## <center> Рекомендательные системы на основе SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данном материале разберёмся с тем, как строить рекомендательные системы на основе сингулярного разложения матриц. \n",
    "\n",
    "Начнём с постановки задачи. Что такое вообще рекомендательная система? Это алгоритм, который предлагает нашим пользователям/клиентам товары/контент с целью увеличения счастья пользователя или выручки.\n",
    "Подходы к построению таких систем могут быть разными:\n",
    "\n",
    "- Рекомендовать самые популярные товары/фильмы/новости/музыку;\n",
    "- Колоборативная фильтрация на основе user-based/item-based (метрические методы);\n",
    "- Рекомендации на основе факторизации матриц(SVD-разложение);\n",
    "- Можно обучить регрессор, предсказывающий выручку или рейтинг фильма для пары \"Пользователь-Объект\";\n",
    "\n",
    "Я же подробно остановлюсь на третьем пункте, т.к. этот метод построения рекомендательных систем достаточно неплох: быстро работает и показывает неплохое качество рекомендаций, выявляя скрытые интересы пользователей и качества объектов (фильм, товар и т.п.)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Представим, что у нас есть небольшой онлайн-кинотеатр и мы хотим, чтобы пользователи были довольны и смотрели фильмы у нас, а не на других ресурсах. Создадим небольшую матрицу рейтингов наших фильмов, на основании оценок пользователей."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Названия фильмов\n",
    "f = [\n",
    "    \"Терминатор\",\n",
    "    \"Робокоп\",\n",
    "    \"Рэмбо\",\n",
    "    \"Джеймс Бонд\",\n",
    "    \"Властелин Колец\",\n",
    "    \"Хоббит\",\n",
    "    \"Гарри Поттер\",\n",
    "    \"Американский Пирог\",\n",
    "    \"Мальчишник в Вегасе\",\n",
    "    \"Мстители\",\n",
    "    \"Супермен\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Наши пользователи\n",
    "u = [\"Вася\", \"Петя\", \"Саша\", \"Женя\", \"Маша\", \"Оля\", \"Лена\", \"Ваня\", \"Ира\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Рейтинг фильмов от пользователей. 0 - значит фильм не был просмотрен.\n",
    "fu = np.array(\n",
    "    [\n",
    "        [9, 0, 0, 6, 4, 3, 5, 7, 8, 6, 5],\n",
    "        [0, 0, 0, 8, 10, 0, 8, 5, 6, 2, 0],\n",
    "        [8, 0, 5, 7, 4, 0, 5, 7, 8, 8, 5],\n",
    "        [5, 2, 0, 4, 10, 0, 9, 5, 0, 0, 0],\n",
    "        [0, 0, 2, 5, 7, 0, 0, 9, 0, 4, 0],\n",
    "        [2, 3, 0, 0, 0, 7, 9, 0, 5, 0, 3],\n",
    "        [5, 0, 3, 8, 8, 0, 0, 7, 0, 2, 4],\n",
    "        [8, 7, 0, 6, 5, 4, 0, 8, 0, 10, 0],\n",
    "        [0, 0, 2, 5, 10, 0, 0, 0, 9, 0, 3],\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Сделаем датафрейм из f,u и fu.\n",
    "R = pd.DataFrame(fu, u, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вот теперь и встаёт вопрос, какой предложить фильм, чтобы наш пользователь хорошо провёл вечер? Немного поискав на просторах интернета, находим статью о том, что один большой онлнайн-кинотеатр **Netflix** в 2006 году объявил соревнование по разработке рекомендательной системы фильмов. За решение, которое улучшит качество рекомендаций на 10%, была объявлена награда в 1 млн долларов. В итоге в 2011 году такое решение было получено, а работало оно на основе **SVD** - Singular Value Decomposition или Сингулярном Разложении Матриц, если по-русски."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Немного линейной алгебры.\n",
    "Смысл SVD состоит в том, что прямоугольную матрицу A мы можем разложить на 3 матрицы: U, Ʃ и V.\n",
    "\n",
    "<img src=\"../../img/SVD_img_2.png\" style=\"width: 300px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Матрицы U и V ортогональные, а Ʃ — диагональная (хотя и не квадратная). \n",
    "\n",
    "<img src=\"../../img/SVD_img_3.png\" style=\"width: 600px\">\n",
    "\n",
    "В случае разреженной матрицы, которой обычно является рейтинг фильмов, можно использовать усечённое разложение (truncated SVD), где мы оставляем только d первых чисел $\\lambda$. В итоге получаем разложение матрицы A', которое хорошо приближает исходную матрицу A."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=../../img/SVD_img_5.png style=\"width: 300px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь применим это к нашей рекомендательной системе. У нас есть матрица оценок R, давайте сделаем её сингулярное разложение.\n",
    "Первые две матрицы перемножим, получим также матрицу $n$ x $d$, в итоге получится такое разложение.\n",
    "<img src=\"../../img/SVD_img_4.png\" style=\"width: 600px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Размерность **d** в данном случае отвечает за количество скрытых свойств/интересов у наших пользователей и фильмов, $\\hat r \\tiny ui$ - рейтинг фильма, который мы получили перемножив матрицы **U** и **V**.\n",
    "\n",
    "Таким образом, мы научились приближать исходную матрицу оценок матрицей, полученной на основе перемножения двух матриц, отвечающих за интересы пользователей и свойства фильмов. Но проблема в том, что матрица **R** у нас разреженная, и мы бы хотели бы как раз заполнить пустые (нулевые) элементы оценками пользователей.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как же это сделать? Для этого нам нужно приблизить уже имеющиеся оценки числами, полученными в результате перемножения матриц $U$ и $V$. Тогда мы получим и оценки для нулей нашей исходной матрицы, т.е. нужно минимизировать такую функцию:\n",
    "\n",
    "$$\\large \\frac{1}{n}\\Sigma_{ui} (r {\\tiny ui}-\\hat r {\\tiny ui})^{2} \\to min,$$ где суммирование идёт по ненулевым индексам $u,i$ нашей исходной матрицы оценок, а $n$ - количество оценок в матрице."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так это же задача регрессии с квадратичной функций потерь **MSE**. **MSE** выбрана по двум причинам: историческая, Netflix предложили метрику **RMSE** и она лучше оптимизируется градиентным спуском."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Функция потерь есть, метод оптимизации выбран, что забыли? Регуляризатор! Воспользуемся $L{\\tiny 2}$-регуляризацией: \n",
    "$\\large \\lambda(\\Sigma_{u}p {\\tiny u}^{2} + \\Sigma_{i}q {\\tiny i}^{2})$\n",
    "Осталось посчитать градиенты для реализации поиска минимума **MSE**. Воспользуемся стохастическим градиентным спуском, поэтому будем считать градиенты на одном объекте $\\large r {\\tiny ui}$ из матрицы оценок. В итоге получаем такие правила обновления элементов наших матриц **U** и **V**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\large \n",
    "{p_{u,j}=p_{u,j} + \\large \\gamma (r_{ui}-\\hat r_{ui}) q_{i,j} - \\lambda p_{u,j}}$\n",
    "\n",
    "\n",
    "$\\large {q_{i,j}=q_{i,j} + \\large \\gamma (r_{ui}-\\hat r_{ui}) p_{u,j} - \\lambda q_{i,j}}$\n",
    "\n",
    "Индексы $u,i$ - номера пользователя и фильма в наших матрицах **U** (номер строки) и **V** (номер столбца), $j$ - $j$-ая компонента векторов $\\large p_{u}$ и $\\large q_{i}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ну что же, теперь воспользуемся Питоном и сделаем рекомендации для наших пользователей!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SVD(R, d, step, lambda_reg, n_iters):\n",
    "    # инициализуем наши матрицы для разложения\n",
    "    U = np.zeros((R.shape[0], d))\n",
    "    V = np.zeros((d, R.shape[1]))\n",
    "    # начальные элементы матрица U и V будут средним рейтингом по ненулевым оценкам\n",
    "    mu = R.sum() / (R != 0).sum()\n",
    "    non_zero = (R != 0).sum()\n",
    "    U = U + mu\n",
    "    V = V + mu\n",
    "    # Создадим списки, где будут индексы нулевых и ненулевых элементов матрицы R\n",
    "    indx = []\n",
    "    zero_indx = []\n",
    "    # Инициализируем MSE в начале и будем отслеживать в процессе обучения\n",
    "    MSE_start = 0\n",
    "    MSE = []\n",
    "    # Найдём индексы нулевых и ненулевых элементов\n",
    "    for i in range(R.shape[0]):\n",
    "        for j in range(R.shape[1]):\n",
    "            if R[i][j] > 0:\n",
    "                indx.append([i, j])\n",
    "                MSE_start += ((R[i, j] - np.dot(U[i, :], V[:, j])) ** 2) / non_zero\n",
    "            else:\n",
    "                zero_indx.append([i, j])\n",
    "    # Сделаем градиентный спуск\n",
    "    for n in range(n_iters):\n",
    "        choice = np.random.randint(0, len(indx))\n",
    "        ij = indx[choice]\n",
    "        for k in range(0, d):\n",
    "            U[ij[0], k] = U[ij[0], k] + step * (\n",
    "                (R[ij[0]][ij[1]] - np.dot(U[ij[0], :], V[:, ij[1]])) * V[k, ij[1]]\n",
    "                - lambda_reg * U[ij[0], k]\n",
    "            )\n",
    "            V[k, ij[1]] = V[k, ij[1]] + step * (\n",
    "                (R[ij[0]][ij[1]] - np.dot(U[ij[0], :], V[:, ij[1]])) * U[ij[0], k]\n",
    "                - lambda_reg * V[k, ij[1]]\n",
    "            )\n",
    "\n",
    "        L = 0\n",
    "        for i in range(R.shape[0]):\n",
    "            for j in range(R.shape[1]):\n",
    "                if R[i, j] > 0:\n",
    "                    L += ((R[i, j] - np.dot(U[i, :], V[:, j])) ** 2) / non_zero\n",
    "        MSE.append(L)\n",
    "    return U, V, MSE_start, MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, B, M_1, M_end = SVD(R.values, 2, 0.01, 0.1, 3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_cap = np.zeros((R.shape[0], R.shape[1]))\n",
    "for i in range(R.shape[0]):\n",
    "    for j in range(R.shape[1]):\n",
    "        r_cap[i, j] = np.dot(A[i, :], B[:, j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_cap = pd.DataFrame(r_cap, u, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_cap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Start MSE:\", M_1, \"Finish MSE:\", M_end[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Start RMSE\", np.sqrt(M_1), \"Finish RMSE\", np.sqrt(M_end[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = (R.values != 0).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 0\n",
    "for i in range(R.shape[0]):\n",
    "    for j in range(R.shape[1]):\n",
    "        if R.values[i, j] > 0:\n",
    "            L += ((R.values[i, j] - R_cap.values[i, j]) ** 2) / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В итоге получили матрицу оценок с помощью SGD. Как видно, MSE сильно упал. Можно поиграться с параметрами, чтобы настроить качество. Как итог, выбираем фильмы с наибольшей оценкой и рекомендуем их пользователям! :) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В заключение хочется сказать, что есть хорошая библиотека Surprise, где можно делать рекомендации проще, быстрее и различными способами.\n",
    "Оставлю ссылку для ознакомления http://surpriselib.com/"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
